Observation of bulk quadrupole in topological heat transport

The quantized bulk quadrupole moment has so far revealed a non-trivial boundary state with lower-dimensional topological edge states and in-gap zero-dimensional corner modes. In contrast to photonic implementations, state-of-the-art strategies for topological thermal metamaterials struggle to achieve such higher-order hierarchical features. This is due to the absence of quantized bulk quadrupole moments in thermal diffusion fundamentally prohibiting possible band topology expansions. Here, we report a recipe for generating quantized bulk quadrupole moments in fluid heat transport and observe the quadrupole topological phases in non-Hermitian thermal systems. Our experiments show that both the real- and imaginary-valued bands exhibit the hierarchical features of bulk, gapped edge and in-gap corner states—in stark contrast to the higher-order states observed only on real-valued bands in classical wave fields. Our findings open up unique possibilities for diffusive metamaterial engineering and establish a playground for multipolar topological physics.


Supplementary Note 1. Effective square lattice for the quadrupole topological phases in heat transport
We first start with the heat exchanges and corresponding heat transport processes of one site described by Eq. (1) of the main content. The right terms of Eq. (1) respectively indicate the conductive and convective processes, as well as the total heat exchanges between the target and its neighbors. Since the heat exchanges are induced by the convective processes, the total exchanged energies can be described by Newton's cooling law, i.e., ∑ / = ∑ ℎ / • ∆ , where h denotes the convective heat transfer coefficient. For the current model, each target site (Ti,j, light-red dashed border) connects four neighbors in the system (Fig. S1a), indicating four independent heat exchange terms related to the target site. Among these heat exchanges, two of them describe the intracell couplings between two neighboring sites within one four-site unit-structure (orange lines of Supplementary Figure 1a), while another two represent the intercell couplings of the heat exchanges between the neighboring sites of two adjacent unit structures (black lines of Supplementary Figure   1a). Due to the different heat exchange capabilities of the intracell and intercell channels, we make ∑ = • ∑ , which can be realized by engineering their convective heat transfer coefficient via modulating the total heat exchange areas of corresponding channels based on Newton's cooling law.
A 2D effective square lattice can be created with the proposed unit-structure, which possesses two advective components respectively along x and y directions. Such a 2D thermal lattice further leads to a thermal grid system (Supplementary Figure 1a), and its energy equation can be expressed as where Ω = Ω / cos( ) and Ω = Ω / sin( ). R(x) and R(y) denote the directional components of one site. i and j depict the spatial information of each site. For simplification, we make the thermal lattice follow the spatial periodicity of 4a along the x and y directions, where a is the distance between the centers of neighboring sites ( , , x y x y x y ij R E c Due to the translation symmetry of the advective configurations, the gradients of decoupled velocity fields along the two vectors (x) and (y) further result in effective potential fields.
. , , Along direction: . In that case, the effective Bloch wave numbers kt,x and kt,y are dependent on the imposed advections, since the periodic wave properties and effective oscillations of the system are determined by their values. Considering the neighboring sites at the interfaces of neighboring lattices (Fig. 1a), the continuous temperature condition leads to the following relations: 1) , and , = • • .

Supplementary Note 2. Theoretical model for thermal quadrupole topological phases
The thermal process of each unit-structure ( .


Then, the effective Hamiltonian of a four-site unit-structure can now be written as The terms (1 + ± , ) and 1 + ± , denote the effective hopping in imaginary parts, which originate from the dissipative thermal couplings between neighboring sites. The terms ±Ω (cos + sin ) and ±Ω (cos + sin ) induced by the imposed advections act as the on-site energies in real parts. The intercell couplings should follow the Bloch theorem, where kt, x and kt, y are the directional components of the effective Bloch wave number kt.
Solving the eigenvalue problem with the wave-like temperatures, the eigenvalues can be obtained can be observed.
Note that both the advections (Ω and Ω ) and the thermal coupling strength ratio ( ) determine the band structure. This suggests the potential of revealing quadrupole topological phases either in an independent real/imaginary band or a complex             Figure 2b) give birth to the under-coupling processes. In that case, we could define the effective positive and negative couplings with these over-coupling and under-coupling behaviors.

Wannier bands and nested Wannier bands for the proposed two strategies
Since the four-site unit structure proposed in the current work satisfies a quadrupole model, the Wannier bands and

Gapped real spectrum induced by Hermitian advections
For the quadrupole topological phases induced by Hermitian advection (Fig. 2 of The changing advection differences ∆Ω lead to the hierarchical properties only in the real-valued bands. The eigenfrequencies sorted in these real-and imaginary-valued bands as the function of ∆Ω are presented in Supplementary   Figures 5a and b. The real band structures of the first Brillouin zones under specific advections are presented in Supplementary Figure   6a. When the imposed advective magnitudes are same and the directions are opposite (Ω = −Ω ), all the real-valued bands degenerate indicating the gapless structure and the transition between nontrivial and trivial thermal states. When we slightly increase Ω to −1.385Ω to maintain an advective difference, two gaps emerge between four bands, and the second and third bands are still degenerate. Further increasing Ω to −3.154Ω , all degeneracies lift and a complete open band emerges. In such transitions, we only need to capture the parity inversion between the Γ and X points of the Brillouin zone due to the C4 symmetry of the square lattice. Here, we consider the lowest band gap (Fig. 1c). Parity inversion at the X point is observed as indicated by the flip of the "+" and "-" signs, thus implying the topological phase transition via solely modulating the Hermitian advection represents a class of topological quadrupole phases, embracing the in-gap 0D and gapped 1D topological modes.

Gapped imaginary spectrum induced by non-Hermitian thermal coupling
For the quadrupole topological phases induced by non-Hermitian thermal coupling (Fig. 4 of The varied coupling strengths of the intercell and intracell channels (β ≠ 1) lead to significant hierarchical properties in imaginary-valued bands. The eigenfrequencies sorted in these real-and imaginary-valued bands as the function of β are illustrated in Supplementary Figures 5c and d

Numerical quadrupole topological phases of the proposed cases
The numerical temperature distributions and experimental field intensities corresponding to Figs. 2 and 4 are illustrated in Supplementary Figures 7 and 8. Significant features of corner, edge, and bulk states are observed in these validations, which overlap well with the experimental findings in the main contents.
It is noting that the observed behaviors are implemented in fluid transport with a non-vanishing first-order term in the energy equation of the N-S equations (advections). In that case, the temperature fields are available to propagate along the advections (drift-terms) which is non-existing in the omni-directional diffusion of pure conduction. Thus, we can easily capture the behaviors on their temperature distributions. The remarkable differences between their behaviors in temperature fields can be reflected on the irregular and regular distributions ( Supplementary Figures 7 and 8), since the convective temperature fields would also carry the information of the inhomogeneous advective vectors. Considering the imposed velocities via the advective balls, the fluid dynamics (transport) within the entire system (inside and outside the advective balls) are also actuated based on the momentum equation + ( • ∇) = −∇ + ∇ + and exhibit these nontrivial behaviors in the real vector space of pressure and velocity. More specifically for the real-valued band realized by the advections in fluid heat transport, the hierarchical states modulated by advective hermiticity in Fig. 2 can be also exhibited with the distributions of pressures or velocity (Supplementary Figure 9).

Robustness of the hierarchical states
The  Figure 11). In that case, the nontrivial states vanish when only modulating these parameters without changing others accordingly. Once modulating these critical values and corresponding parameters simultaneously to satisfy the frequency shifts, the nontrivial states reoccur again and maintain robustly.
Due to the non-zero nested Wannier bands for these two strategies, the robustness of these hierarchical states is also independent of the site numbers when the critical parameters for creating the bands are unchanged. Additional validations respectively with 192 and 96 sites are implemented (Supplementary Figure 12), whose input parameters are same with the cases in Figs. 2 and 4 of the main content. Significant corner states are observed with the two strategies regardless of the site numbers.

Robustness of the hierarchical states under defects
We further validate the robustness of these hierarchical states under defects of the lattice. To create the defects, we

Supplementary Note 5. Time changing rate of the temperature intensity
The observed behaviors in the main constant are measured at steady states, after the temperature field evolutions reaching equilibrium. Before reaching the observed states, the systems are nonequilibrium. We found that these hierarchical states are also available under time-dependent dynamics, since the temperature field intensities of corresponding sites are always higher than the others at specific times. To characterize such behaviors at nonequilibrium states, we adopt the time changing rate of the field intensity measured by experiments ( , =  Figure 15). The smaller values of imply the small changes in temperature profiles during the nonequilibrium processes, thus revealing the higher intensities at specific sites with energy localizations.

Supplementary Note 6. Potential impacts on general diffusive transports
The realizations of such quadrupole topology in non-Hermitian heat transport indicate a distinct route to modulate thermal behaviors and maintain robustness under changing and complex ambient. Rapidly increasing energy density and heat flux within the thermal techniques in almost all industrial fields, such as conventional energy power and emerging electronics, raise new demands of higher-efficiency thermal energy utilization, and more intelligent thermal management.
Both two aspects directly touch the most fundamental diffusive nature of heat transport, which notoriously restricts the dynamic coherences in principle, let alone the robust heat control and multifarious field reconfigurations. The different and robust profiles of heat transports observed in the work might be a perfect strategy for solving these challenging problems and further pave the way of flexible manipulations. Based on the observed heat patterns, the most intuitive application is heat management, since we can easily localize heat energies at tailored positions of the system independent of the materials and structures. Further exploring the merging states, the localized positions can be anywhere in the system, which is not limited to the bulk, edge, and corners. Such behaviors can be generally realized by creating an internal bound state within the system 1 , which break the lattice translational symmetry. As an example, some state-of-art works exploit the butting of topological nontrivial and trivial lattices to form "long-long defects" at their boundaries 1 . Such a method requires the combinations of nontrivial and trivial lattice within the entire system, thus naturally break the translational symmetry at their boundaries. Another potential approach is to insert domain walls at tailored sites (Supplementary Figure 16) in the system possessing one lattice type (such as the square-lattice of the current work). In that case, the neighboring lattices at the interface are mismatched thus giving rise to the internal bound states. Such a property deserves more studies, and it is quite important for exploiting paradigm-shift thermal management in electronic devices and other industrial applications.
It is worth noting that heat transport is one of the typical information carried by fluid dynamics and the extensive transport phenomena, thus we can also reveal similar behaviors in other transport fields under the same framework proposed in this work, and further shed new lights on the unexpected control of mass, charge, momentum, and fluid fields.
Such universality of the proposed theory is enabled by the same form of constitutive equations for different transports, such as the Fick low for mass and particle transports, the viscosity for momentum conversion, the drift-diffusion equation